home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / MATH1 / ERFD3.PAS < prev    next >
Pascal/Delphi Source File  |  1985-04-03  |  1KB  |  83 lines
  1. program erfd3;        { -> 330 }
  2. { evaluation of the gaussian error function }
  3.  
  4. var    x,er,ec        : real;
  5.     done        : boolean;
  6.  
  7. external procedure cls;
  8.  
  9. function erf(x: real): real;
  10. { infinite series expansion of the Gaussian error function }
  11.  
  12. const    sqrtpi        = 1.7724538;
  13.     tol        = 1.0E-4;
  14.  
  15. var    x2,sum,sum1,term: real;
  16.     i        : integer;
  17.  
  18. begin
  19.       x2:=x*x;
  20.       sum:=x;
  21.       term:=x;
  22.       i:=0;
  23.       repeat
  24.     i:=i+1;
  25.     sum1:=sum;
  26.     term:=2.0*term*x2/(1.0+2.0*i);
  27.     sum:=term+sum1
  28.       until term<tol*sum;
  29.       erf:=2.0*sum*exp(-x2)/sqrtpi
  30. end;    { erf }
  31.  
  32. function erfc(x: real): real;
  33. { complement of error function }
  34. const    sqrtpi        = 1.7724538;
  35.     terms        = 12;
  36.  
  37. var    x2,u,v,sum    : real;
  38.     i        : integer;
  39. begin
  40.   x2:=x*x;
  41.   v:=1.0/(2.0*x2);
  42.   u:=1.0+v*(terms+1.0);
  43.   for i:=terms downto 1 do
  44.     begin
  45.       sum:=1.0+i*v/u;
  46.       u:=sum
  47.     end;
  48.   erfc:=exp(-x2)/(x*sum*sqrtpi)
  49. end;        { ercf }
  50.  
  51. begin        { main }
  52.   cls;
  53.   done:=false;
  54.   writeln;
  55.   repeat
  56.     write('Arg? ');
  57.     readln(x);
  58.     if x<0.0 then done:=true
  59.     else
  60.       begin
  61.     if x=0.0 then
  62.       begin
  63.         er:=0.0;
  64.         ec:=1.0
  65.       end
  66.     else
  67.       begin
  68.         if x<1.5 then
  69.           begin
  70.         er:=erf(x);
  71.         ec:=1.0-er
  72.           end
  73.         else
  74.           begin
  75.         ec:=erfc(x);
  76.         er:=1.0-ec
  77.         end    { if }
  78.     end;
  79.     writeln('X= ',x,' Erf= ',er:7:4,', Erfc= ',ec)
  80.       end    { if }
  81.     until done
  82. end.
  83.